Actividad unidad 2

Problema 1

Estimación del valor de \(\pi\)

La siguiente figura sugiere como estimar el valor de \(\pi\) con una simulación. En la figura, un circuito con un área igual a \(\pi/4\), está inscrito en un cuadrado cuya área es igual a 1. Se elige de forma aleatoria \(n\) puntos dentro del cuadrado. La probabilidad de que un punto esté dentro del círculo es igual a la fracción del área del cuadrado que abarca a este, la cual es \(\pi/4\). Por tanto, se puede estimar el valor de \(\pi/4\) al contar el número de puntos dentro del círculo, para obtener la estimación de \(\pi/4\). De este último resultado se encontrar una aproximación para el valor de \(\pi\).

Imagen con circulo circunscrito

Pasos sugeridos:

Genere n coordenadas \(x: X1, . . . , Xn\). Utilice la distribución uniforme con valor mínimo de 0 y valor máximo de 1. La distribución uniforme genera variables aleatorias que tienen la misma probabilidad de venir de cualquier parte del intervalo (0,1).

Genere 1000 coordenadas \(y: Y1,...,Yn\), utilizando nuevamente la distribución uniforme con valor mínimo de 0 y valor máximo de 1.

Datos simulación

set.seed(123)
n1 <- 1000
n2 <- 10000
n3 <- 100000
x_1000 <- runif(n1, 0, 1)
y_1000 <- runif(n1, 0, 1)

x_10000 <- runif(n2, 0, 1)
y_10000 <- runif(n2, 0, 1)

x_100000 <- runif(n3, 0, 1)
y_100000 <- runif(n3, 0, 1)

Cada punto (Xi,Yi) se encuentra dentro del círculo si su distancia desde el centro (0.5,0.5) es menor a 0.5. Para cada par (Xi,Yi) determine si la distancia desde el centro es menor a 0.5. Esto último se puede realizar al calcular el valor (Xi−0.5)2+(Yi−0.5)2, que es el cuadrado de la distancia, y al determinar si es menor que 0.25.

Respuestas

¿Cuántos de los puntos están dentro del círculo?

puntos_dentro_1000 <- sum((x_1000 - 0.5)^2 + (y_1000 - 0.5)^2 <= 0.25)
puntos_dentro_10000 <- sum((x_10000 - 0.5)^2 + (y_10000 - 0.5)^2 <= 0.25)
puntos_dentro_100000 <- sum((x_100000 - 0.5)^2 + (y_100000 - 0.5)^2 <= 0.25)

El número de puntos dentro del círculo es de: 800

¿Cuál es su estimación de \(\pi\)?

\(\pi = 4 \times \frac{\text{puntos dentro}}{n}\)

pi_estimado_1000 <- 4 * puntos_dentro_1000 / n1
pi_estimado_10000 <- 4 * puntos_dentro_10000 / n2
pi_estimado_100000 <- 4 * puntos_dentro_100000 / n3

La estimación de \(\large\pi\) con \(\large n=1000\) es de: 3.2

La estimación de \(\large\pi\) con \(\large n=10^{4}\) es de: 3.1528

La estimación de \(\large\pi\) con \(\large n=10^{5}\) es de: 3.144

¿Cuál es el error absoluto de su estimación?

\(\Large|\pi - \pi_{estimado}|\)

error_absoluto_1000 <- abs(pi - pi_estimado_1000)
error_absoluto_10000 <- abs(pi - pi_estimado_10000)
error_absoluto_100000 <- abs(pi - pi_estimado_100000)

Error absoluto para muestra de \(n=1000\) : \(0.0584073\)

Error absoluto para muestra de \(n=10^{5}\) : \(0.0112073\)

Error absoluto para muestra de \(n=10^{4}\) : \(0.0024073\)

¿Cuál es el error relativo de su estimación?

\(\Large\frac{|\pi - \pi_{estimado}|}{\pi} = \frac{|\pi - 3.2|}{\pi}\)

error_relativo_1000 <- abs(pi - pi_estimado_1000) / pi
error_relativo_10000 <- abs(pi - pi_estimado_10000) / pi
error_relativo_100000 <- abs(pi - pi_estimado_100000) / pi

Error relativo para muestra de \(n=1000\) : \(0.0185916\)

Error relativo para muestra de \(n=10^{5}\) : \(0.0035674\)

Error relativo para muestra de \(n=10^{4}\) : \(7.6628216\times 10^{-4}\)

Problema 2

Propiedades de los estimadores

La simulación ayuda a entender y validad las propiedades de los estimadores estadísticos como son. insesgadez, eficiencia y la consistencia principalmente. El siguiente problema permite evidenciar las principales características de un grupo de estimadores propuestos para la estimación de un parámetro asociado a un modelo de probabilidad.

Sean \(X1, X2, X3 \,y\, X4\), una muestra aleatoria de tamaño \(n=4\) cuya población la conforma una distribución exponencial con parámetro \(\theta\) desconocido. Determine las características de cada uno de los siguientes estimadores propuestos:

\(\hat{\theta}_1 = \frac{X_1 + X_2}{6} + \frac{X_3 + X_4}{3}\)
\(\hat{\theta}_2 = \frac{X_1 + 2X_2 + 3X_3 + 4X_4}{5}\)
\(\hat{\theta}_3 = \frac{X_1 + X_2 + X_3 + X_4}{4}\)
\(\hat{\theta}_4 = \frac{\min\{X_1, X_2, X_3, X_4\} + \max\{X_1, X_2, X_3, X_4\}}{2}\)

Formulación de la simulación

n <- 4
m = 5000
lam <- 0.5

rep = rexp(m*n,rate = lam)
df = data.frame(matrix(rep,nrow = m, ncol = n, byrow = TRUE))

estimador1 <- function(m) { ((m[1] + m[2]) / 6) + ((m[3] + m[4]) / 3) }
estimador2 <- function(m) { (m[1] + 2*m[2] + 3*m[3] + 4*m[4]) / 5 }
estimador3 <- function(m) { mean(m) }
estimador4 <- function(m) { (min(m) + max(m)) / 2 }

estimarMuestra <- function(data, n_muestra, lam) {
  teta <- 1 / lam
  muestra <- data[sample(nrow(df), size=n_muestra), ]
  t1 <- apply(muestra, 1, estimador1)
  t2 <- apply(muestra, 1, estimador2)
  t3 <- apply(muestra, 1, estimador3)
  t4 <- apply(muestra, 1, estimador4)
  return(data.frame(t1, t2, t3, t4))
}
calcularMetricas <- function(datos, n_muestra, lam) {
  teta <- 1 / lam
  media <- apply(datos, 2, mean)
  varianza <- apply(datos, 2, var)
  sesgo <- media - teta
  return (data.frame(media = media, varianza = varianza, n=n_muestra, sesgo=sesgo))
}

obtenerGrafica <- function(datos, n_muestra, lam) {
  teta <- 1 / lam
  d <- pivot_longer(datos, cols = c(t1, t2, t3, t4), names_to = "categoria", values_to = "valor")
  bp <- ggplot(d, aes(x = categoria, y = valor, fill = categoria)) +
    geom_boxplot() +
    geom_hline(yintercept = teta, color = "red", linetype = "dashed", linewidth = 1) +
    scale_fill_manual(values = c("#377eb8", "#ff7f00", "#4daf4a", "#f781bf")) +
    labs(title = paste("m =", n_muestra), x = "Categoría", y = "Valor") +
    theme_minimal() 
  bp <- ggplotly(bp)
  return (bp)
}

Resultados m=20

Resultados m=50

Resultados m=100

Resultados m=1000

Conclusiones

Podemos ver que \(\hat\theta_1\) y \(\hat\theta_3\) son insesgados, eficientes y consistentes para las muestras \(n = 20\), \(n = 50\),\(n = 100\) y \(n = 1000\).
Para \(\hat\theta_4\) es insesgado, eficiente y consistente para las muestras \(n > 100\).
El estimador \(\hat\theta_2\) para sus tres indicadores es muy grande, aunque se ve que comienza a mejorar a medida que el tamaño de la muestra aumenta.

Problemas 3

Teorema del Límite Central

El Teorema del Límite Central es uno de los más importantes en la inferencia estadística y habla sobre la convergencia de los estimadores como la proporción muestral a la distribución normal. Algunos autores afirman que esta aproximación es bastante buena a partir del umbral \(n>30\).

A continuación se describen los siguientes pasos para su verificación:

Realice una simulación en la cual genere una población de \(n=1000\) (Lote), donde el porcentaje de individuos (supongamos plantas) enfermas sea del 50%.
Genere una función que permita obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral \(\hat p\) para un tamaño de muestra dado \(n\).
Repita el escenario anterior \((b)n=500\) veces y analice los resultados en cuanto al comportamiento de los 500 resultados del estimador \(\hat p\). ¿Qué tan simétricos o sesgados son los resultados obtenidos? y ¿qué se puede observar en cuanto a la variabilidad?. Realice en su informe un comentario sobre los resultados obtenidos.
Repita los puntos b y c para tamaños de muestra \(n=5, 10, 15, 20, 30, 50, 60, 100, 200, 500\). Compare los resultados obtenidos para los diferentes tamaños de muestra en cuanto a la normalidad. Utilice pruebas de bondad y ajuste (shapiro wilks :shspiro.test()) y métodos gráficos (gráfico de normalidad: qqnorm()). Comente en su informe los resultados obtenidos.
Repita toda la simulación (puntos a – d), pero ahora para lotes con 10% de plantas enfermas y de nuevo para lotes con un 90% de plantas enfermas. Concluya sobre los resultados del ejercicio.

Formulación de la simulación

a. Generando muestra 50-50

generador_valores <- function(n, porcentaje1, procentaje2) {
  p_sano <- round(n * porcentaje1)
  p_enfermo <- round(n * procentaje2)
  estado <- c(rep("sana", p_sano), rep("enferma", p_enfermo))
  return(estado)
}
lote <- generador_valores(1000, 0.5, 0.5)

b. funcion que permita obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral \(\hat p\) para un tamaño de muestra dado \(n\).

generarMuestra <- function(lote, n) {
  muestra <- sample(lote, n)
  return(muestra)
}
calcularProporcion <- function(lote, n, estado) {
  muestra <- generarMuestra(lote, n)
  P <- sum(muestra == estado) / n
  return(P)
}

proporcion <- calcularProporcion(lote, 10, "enferma")

set.seed(123)
n <- 500
resultados <- replicate(500, calcularProporcion(lote, 100, "enferma"))
hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")

¿Qué tan simétricos o sesgados son los resultados obtenidos?

¿qué se puede observar en cuanto a la variabilidad?

Repita los puntos b y c para

iterador <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
set.seed(123)
for(n in iterador) {
  resultados <- replicate(n, calcularProporcion(lote, 100, "enferma"))
  par(mfrow=c(1,2))
  hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
  qqnorm(resultados, col = "#377eb8")
  qqline(resultados, col = "#ff7f00")
}

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.89512, p-value = 0.3835

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.88478, p-value = 0.148

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.93633, p-value = 0.3385

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.95858, p-value = 0.5158

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98693, p-value = 0.9654

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.95157, p-value = 0.03958

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97943, p-value = 0.405

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.99032, p-value = 0.6906

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98887, p-value = 0.1217

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.9949, p-value = 0.09726

set.seed(123)
lote <- generador_valores(1000, 0.9, 0.1)
iterador <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
for(n in iterador) {
  resultados <- replicate(n, calcularProporcion(lote, 100, "enferma"))
  par(mfrow=c(1,2))
  hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  qqnorm(resultados, col = "#377eb8", main = paste("QQ de la muestra n =", n))
  qqline(resultados, col = "#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
}

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.95235, p-value = 0.754

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.87409, p-value = 0.1115

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.9518, p-value = 0.5532

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97029, p-value = 0.761

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.94864, p-value = 0.1555

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97431, p-value = 0.3435

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97956, p-value = 0.4102

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.9808, p-value = 0.1535

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98369, p-value = 0.02036

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.9844, p-value = 3.366e-05

Repita toda la simulación (puntos a – d),

set.seed(123)
lote <- generador_valores(1000, 0.1, 0.9)
iterador <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)

for(n in iterador) {
  resultados <- replicate(n, calcularProporcion(lote, 100, "enferma"))
  par(mfrow=c(1,2))
  hist(resultados, col="#377eb8", main = paste("Histograma de la muestra n =", n), xlab="Probabilidad de las muestras", ylab="Frecuencia", las=1, font.axis=4)
  abline (v=mean(resultados), lwd = 4, lty = 2, col="#ff7f00")
  shapiroRes <- shapiro.test(resultados)
  print(shapiroRes)
  qqnorm(resultados, col = "#377eb8")
  qqline(resultados, col = "#ff7f00")
}

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.78818, p-value = 0.06469

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97435, p-value = 0.928

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.94281, p-value = 0.4191

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.89915, p-value = 0.03975

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.94828, p-value = 0.152

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97632, p-value = 0.4092

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.96018, p-value = 0.04802

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.97212, p-value = 0.03221

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.96782, p-value = 0.0001535

## 
##  Shapiro-Wilk normality test
## 
## data:  resultados
## W = 0.98223, p-value = 8.595e-06

Actividad unidad 2

2024-03-02

Problema 1

Estimación del valor de \(\pi\)

Datos simulación

Respuestas

Problema 2

Propiedades de los estimadores

Formulación de la simulación

Resultados m=20

Resultados m=50

Resultados m=100

Resultados m=1000

Conclusiones

Problemas 3

Teorema del Límite Central

Formulación de la simulación

a. Generando muestra 50-50

b. funcion que permita obtener una muestra aleatoria de la población y Calcule el estimador de la proporción muestral \(\hat p\) para un tamaño de muestra dado \(n\).

Repita los puntos b y c para

Repita toda la simulación (puntos a – d),